% =====================================================================
% ====  "Risk-Sharing in Village Economies Revisited", by Bold and Broer

% ====   This programme choses the preference parameters (CRRA rho,
% discount factor beta) according to various criteria

% ====    Written by Tessa Bold and Tobias Broer

% ====    This version June 2021
% =====================================================================%
clear
clf


      
Momvec1data=NaN(18,1);
Momvec2data=NaN(18,1);
Tessa1Tobi2=1;
outlier=0;
analytic=0;
est_rho1=0;
unconditional=0;
village=1
relvarincscale=1;
COALDEV=[3];
coaldev=3;
% =====================================================================%
% ====    Data moments
% =====================================================================%

clear icrisatdata momentsdatabootstrap DCdata DYdata DCdataneg DYdatapos DCdataneg DYdataneg Cdata Ydata VCV1data VCV2data
clear RelCYVardata RelDCDYVardata betaDCDYdatadyneg betaDCDYdatadypos betaDCDYdata Momvec1datamat Momvec2datamat DCskewdata RelDCDYVardata Cdata Ydata DCdata DYdata
clear betaDCDYdata CVARdata YVARdata DCVARdata DYVARdata Cskewdata DCskewdata RelCYVardata RelDCDYVardata


    
if Tessa1Tobi2==1
    outputpath='C:\Users\tbold\Dropbox\Tessa and Tobi\JEEA Replication\Programmes\Main Results and Robustness';
    datapath='C:\Users\tbold\Dropbox\Tessa and Tobi\Data/Output';
    programpath='C:\Users\tbold\Dropbox\Tessa and Tobi\JEEA Replication\Programmes';
else
  
    outputpath='C:\Users\tbroer\Dropbox\Research\Current_Projects\tessa\JEEA Replication\Programmes\Main Results and Robustness';
    datapath='C:\Users\tbroer\Dropbox\Research\Current_Projects\tessa\Data/Output';
    programpath='C:\Users\tbroer\Dropbox\Research\Current_Projects\tessa\JEEA Replication\Programmes';
end



    datafile=['icrisatmoments' int2str(outlier) ];
 
    cd(datapath)
    
    eval(['load ' datafile '.out'])
    eval(['datamoments=' datafile ';']);
    eval(['clear ' datafile])
    
    
    Nom=4;
    
    bootstrapvar=datamoments(1:Nom,4+(village-1)*Nom+1:4+(village)*Nom);
    analyticvar=datamoments(1:Nom,16+(village-1)*Nom+1:16+(village)*Nom);
    analyticmoments=datamoments(1:Nom,1:3);
    
    analyticrawmoments=datamoments(1:Nom,28+1:28+3);
    bootstraprawvar=datamoments(1:Nom,28+4+(village-1)*Nom+1:28+4+(village)*Nom);
    analyticrawvar=datamoments(1:Nom,28+16+(village-1)*Nom+1:28+16+(village)*Nom);
    %%%here you replace the exact standard errors (the variances are
    %%%bootstrapped in STATA, the coefficients are delta-method
    
    VCV2conddata=bootstrapvar;
    VCV2rawdata=bootstraprawvar;
    
      
    
    Momvec2conddata=analyticmoments(:,village)';
    Momvec2rawdata=analyticrawmoments(:,village)';
    if unconditional==1
        VCV2data=VCV2rawdata;
        Momvec2data=Momvec2rawdata;
    else
        VCV2data=VCV2conddata;
        Momvec2data=Momvec2conddata;  
    end
    
    COUNT=0;

        
        cd(outputpath)
            
        eval(['load ROT_Momentsest.mat'])
        
        
        
        
        betavec=[0.94:0.005:0.99];
  
        
        
        
    
               
            
       for inddelta=1:length(betavec)
                for soi=1:6
                    for quasi=1:3
                    rr=find(isfinite(Momentsest(1,1,inddelta,1,:,soi,quasi)) & Momentsest(1,1,inddelta,1,:,soi,quasi)>0 ) 
                    vss_max=max(rr);
                    Vardylog=Momentsest(22,1,inddelta,1,vss_max,soi,quasi);              
                    Vardylog_cond=Momentsest(30,1,inddelta,1,vss_max,soi,quasi);
                             

                            Momentstable(1:3,inddelta,soi,quasi)=Momentsest(1:3,1,inddelta,1,vss_max,soi,quasi);
                            Momentstable(4:7,inddelta,soi,quasi)=[Momentsest([8 15],1,inddelta,1,vss_max,soi,quasi) ; Momentsest(9,1,inddelta,1,vss_max,soi,quasi)/Vardylog_cond(village) ;(Momentsest([16],1,inddelta,1,vss_max,soi,quasi)-Momentsest([17],1,inddelta,1,vss_max,soi,quasi))];
                            Momentstable(8,inddelta,soi,quasi)=(Momentsest(21,1,inddelta,1,vss_max,soi,quasi)-Momentsest(21,1,inddelta,1,vss_max,soi,quasi))./Momentsest(21,1,inddelta,1,vss_max,soi,quasi);%2*maxcerr*(j-1)/(size(moments,2)-1)/moments(28,j);
                            Momentstable(9,inddelta,soi,quasi)=2*Momentsest(5,1,inddelta,1,vss_max,soi,quasi)./Vardylog;%Momentsest(22,j,inddelta,indrho,vss);
                                             
                            Momentstable(1,inddelta,soi,quasi)=Momentstable(1,inddelta,soi,quasi)+1;
                                                      
                            MOMENTSDIAGNOSTICS(1:9,inddelta,soi,quasi)=Momentstable(1:9,inddelta,soi,quasi);
                            
                            momentsdev=Momentstable(4:7,inddelta,soi,quasi)- Momvec2data(1,1:4)';
                            %%%2,4, d:diagonal, 
                            
                            
                            
                            W4(inddelta,soi,quasi)=momentsdev'*inv(VCV2data(1:4,1:4))*momentsdev;
                            W4d(inddelta,soi,quasi)=momentsdev'*inv(diag(diag(VCV2data(1:4,1:4))))*momentsdev;
                            W4_eye(inddelta,soi,quasi)=momentsdev'*inv(eye(4))*momentsdev;

                            W2(inddelta,soi,quasi)=momentsdev(1:2)'*inv(VCV2data(1:2,1:2))*momentsdev(1:2);
                            W2d(inddelta,soi,quasi)=momentsdev(1:2)'*inv(diag(diag(VCV2data(1:2,1:2))))*momentsdev(1:2);
                                                   
                            
                    end
                end
       end
    %%Estimation  
    for soi=1:6
        for quasi=1:3
     WW=W4d(:,soi,quasi);
     [CC,II]=min(WW);
     solution_4(:,soi,quasi)=[Momentstable(1:9,II,soi,quasi); CC];
        end
    end
    
 %Stationary and full risk-sharing gives the best fit%%%
 solution_4(:,:,1)
 solution_4(:,:,2)
 solution_4(:,:,3)
 
[CC,II]=min(solution_4(end,:,1))
[CC,II]=min(solution_4(end,:,2))
[CC,II]=min(solution_4(end,:,3))





